

#####motifbreakR
library(motifbreakR)
library(BSgenome.Hsapiens.UCSC.hg19)
library(SNPlocs.Hsapiens.dbSNP142.GRCh37)
library(SNPlocs.Hsapiens.dbSNP144.GRCh37)
library(BiocParallel)

file=read.table("bias_AShM_anno.hg19_multianno.csv",head=T,sep=",")
snps=as.character(file$avsnp150)
snps=snps[!snps=="."]

search.genome=BSgenome.Hsapiens.UCSC.hg19
dbSNP=SNPlocs.Hsapiens.dbSNP142.GRCh37
rsid=snps

change.to.search.genome <- function(granges.object, search.genome) {
  if (Reduce("&", !is.na(genome(granges.object)))) {
    if (identical(genome(granges.object), genome(search.genome))) {
      return(granges.object)
    }
  }
  if(isTRUE(all.equal(seqlevels(granges.object), seqlevels(search.genome)))) {
    seqinfo(granges.object) <- seqinfo(search.genome)
  } else {
    if(seqlevelsStyle(granges.object) != seqlevelsStyle(search.genome)) {
      seqlevelsStyle(granges.object) <- seqlevelsStyle(search.genome)
    }
    normal.xome <- seqlevels(granges.object)[(regexpr("_", seqlevels(granges.object)) < 0)]
    #xome.value <- str_extract(normal.xome, "[0-9]|1[0-9]|2[0-9]|3[0-9]|4[0-9]|5[0-9]|6[0-9]|7[0-9]|8[0-9]|9[0-9]|X|Y|M")
    positions <- unlist(sapply(paste0(normal.xome, "$"), grep, seqnames(seqinfo(search.genome))))
    new2oldmap <- rep(NA, length(seqinfo(search.genome)))
    new2oldmap[positions] <- 1:length(positions)
    seqinfo(granges.object, new2old = new2oldmap) <- seqinfo(search.genome)
  }
  return(granges.object)
}

rsid.grange <- as(snpsById(dbSNP, rsid, ifnotfound = "warning"),"GRanges")
rsid.grange <- change.to.search.genome(rsid.grange, search.genome)
rsid.grange <- GRanges(rsid.grange)
rsid.refseq <- getSeq(search.genome, rsid.grange)
rsid.grange$UCSC.reference <- as.character(rsid.refseq)
rsid.grange <- split(rsid.grange, rsid.grange$RefSNP_id)

determine.allele.from.ambiguous <- function(ambiguous.allele, known.allele) {
  neucleotide.ambiguity.code <- list(Y = c("C", "T"), R = c("A", "G"), W = c("A", "T"),
                                     S = c("G", "C"), K = c("T", "G"), M = c("C", "A"),
                                     D = c("A", "G", "T"), V = c("A", "C", "G"),
                                     H = c("A", "C", "T"), B = c("C", "G", "T"),
                                     N = c("A", "C", "G", "T"))
  specnac <- neucleotide.ambiguity.code[[ambiguous.allele]]
  unknown.allele <- specnac[-grep(known.allele, specnac)]
  return(unknown.allele)
}
snps=rsid.grange
rt.alt.allele=c()
rsid2=c()
for (i in rsid){
snp=snps[[i]]
if (!is.null(snp)){
alt.allele <- determine.allele.from.ambiguous(snp$alleles_as_ambig, snp$UCSC.reference)
if (length(alt.allele)>=1){
rt.alt.allele=c(rt.alt.allele,alt.allele)
rsid2=c(rsid2,i)
}
}
}
length(rsid2)

snps=rsid2
snps.mb <- snps.from.rsid(snps,dbSNP = SNPlocs.Hsapiens.dbSNP142.GRCh37,search.genome = BSgenome.Hsapiens.UCSC.hg19)
pwmlist=query (query (MotifDb, 'hsapiens'), 'SELEX')
rt=motifbreakR(snps.mb, pwmList=pwmlist, threshold = 1e-4, filterp = TRUE,method = "log", show.neutral = FALSE, verbose = FALSE,bkg = c(A = 0.25, C = 0.25, G = 0.25, T = 0.25), BPPARAM = bpparam())
a=data.frame(rt)
a$rsid=names(rt)
write.table(a,"lncrnas_merge_regulated_SNPs_motif.txt",quote=F,row.names=F,sep="\t")

###
snps.mb=snps.from.file(file = file,
                                  dbSNP = SNPlocs.Hsapiens.dbSNP142.GRCh37,
                                  search.genome = BSgenome.Hsapiens.UCSC.hg19,
                                  format = "bed")
#####
file1=read.table("1202ASMs_motif.txt",head=T,sep="\t")
file2=read.table("1202ASMs_con_motif2.txt",head=T,sep="\t")
file1=file1[file1$dataSource=="jaspar2016",]
file2=file2[file2$dataSource=="jaspar2016",]
table(file1$effect=="strong")
#file11=file1[file1$effect=="strong",]
#file22=file2[file2$effect=="strong",]
a=data.frame(rsid_case=file1$rsid,geneSymbol=file1$geneSymbol,scoreRef_case=file1$scoreRef,scoreAlt_case=file1$scoreAlt,effect_case=file1$effect)
b=data.frame(rsid_con=file2$rsid,geneSymbol=file2$geneSymbol,scoreRef_con=file2$scoreRef,scoreAlt_con=file2$scoreAlt,effect_con=file2$effect)
c=merge(a,b,by="geneSymbol")
symbol=as.character(unique(c$geneSymbol))
result=data.frame(matrix(NA,length(symbol),ncol=6),row.names = symbol)
names(result)=c("case_W","case_S","con_W","con_S","OR","p.value")
for (i in 1:length(symbol)){
test=c[c$geneSymbol==symbol[i],]
test_a=test[!duplicated(test$rsid_case),]
result[i,1]=unlist(dim(test_a[test_a$effect_case=="weak",]))[1]
result[i,2]=unlist(dim(test_a[test_a$effect_case=="strong",]))[1]
test_b=test[!duplicated(test$rsid_con),]
result[i,3]=unlist(dim(test_b[test_b$effect_con=="weak",]))[1]
result[i,4]=unlist(dim(test_b[test_b$effect_con=="strong",]))[1]
result[i,5]=fisher.test(matrix(c(result[i,1],result[i,2],result[i,3],result[i,4]),nrow=2))$estimate
result[i,6]=fisher.test(matrix(c(result[i,1],result[i,2],result[i,3],result[i,4]),nrow=2))$p.value
}
result$geneSymbol=symbol


